﻿#include <math.h>
#include <iostream>
#include <fstream>
#include <cassert>
#include <sstream>
#include <cstdlib>
#include "gdal.h"
#include "gdal_priv.h"
#include "LandSurfaceTemp.h"

void Landsat8_LST::radImg(float* band, float add, float bias, int imglength) {
	for (int i = 0; i < imglength; i++) {
		band[i] = band[i] * bias - add;
	}
}

float* Landsat8_LST::SurEmiss(float* nir, float* red, int imglength) {
	float* aima = new float[imglength];
	memset(aima, 0, imglength);
	for (int i = 0; i < imglength; i++)
	{
		float NDVI = (nir[i] - red[i]) / (nir[i] + red[i]);
		float FV = 0;
		if (NDVI > 0.7)
		{
			FV = 1;
		}
		if (NDVI < 0)
		{
			FV = 0;
		}
		if ((NDVI < 0.7) && (NDVI > 0))
		{
			FV = (NDVI - 0.05) / (0.7 - 0.05);
		}
		aima[i] = FV * 0.004 + 0.986;
	}

	return aima;
}


float* Landsat8_LST::LST(float t, float Lu, float Ld, float* aima, float* band10_rad, int imglength) {
	float* LST = new float[imglength];
	float* BlaRad = new float[imglength];

	for (int i = 0; i < imglength; i++)
	{
		BlaRad[i] = (band10_rad[i] - Lu - t * (1 - aima[i]) * Ld) / (t * aima[i]);
		LST[i] = 1321.08 / log(774.89 / BlaRad[i] + 1) - 273;
	}
	delete[] BlaRad;
	delete[] aima;
	return LST;
}

void Landsat8_LST::readImage(char* imgpath, ImgData& IMG) {

	GDALDataset* img = (GDALDataset*)GDALOpen(imgpath, GA_ReadOnly);
	if (img != NULL) {

		int imgWidth = img->GetRasterXSize(); 
		int imgHeight = img->GetRasterYSize();
		int bandNum = img->GetRasterCount();

		IMG.imgBandNum = bandNum;
		IMG.imgHeight = imgHeight;
		IMG.imgWidth = imgWidth;

		GDALRasterBand* poBand = img->GetRasterBand(1);
		img->GetGeoTransform(IMG.adfGeoTransform);
		IMG.proj = string(img->GetProjectionRef());

		int size = imgWidth * imgHeight * bandNum;
		IMG.pData = new float[size];
		IMG.imglength = size;

		poBand->RasterIO(GF_Read, 0, 0, imgWidth, imgHeight, IMG.pData,
			imgWidth, imgHeight, GDT_Float32, 0, 0);

		GDALClose(img);
	}
}

void Landsat8_LST::writeImage(char* imgPath, float* Img, ImgData& tpl) {
	GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
	GDALDataset* poDataset = poDriver->Create(imgPath, tpl.imgWidth, tpl.imgHeight, tpl.imgBandNum, GDT_Float32, NULL);
	poDataset->SetGeoTransform(tpl.adfGeoTransform);
	poDataset->SetProjection(tpl.proj.c_str());
	poDataset->RasterIO(GF_Write, 0, 0, tpl.imgWidth, tpl.imgHeight,
		Img, tpl.imgWidth, tpl.imgHeight, GDT_Float32, tpl.imgBandNum, 0, 0, 0, 0);
	GDALClose(poDataset);
}


string Landsat8_LST::removeNUllString(string sNewTag) {
	int begin = 0;
	begin = sNewTag.find(" ", begin);
	while (begin != -1)
	{
		sNewTag.replace(begin, 1, "");
		begin = sNewTag.find(" ", begin);
	}
	return sNewTag;
}

template <class Type>
Type Landsat8_LST::stringToNum(const string & str)
{
	istringstream iss(str);
	Type num;
	iss >> num;
	return num;
}

void Landsat8_LST::getMetaData(string metaDataPath, float* R_add, float* NIR_add, float* TIRS_add,
	float* R_mult, float* NIR_mult, float* TIRS_mult, string getTime,
	float* longtitude, float* latitude) {

	std::ifstream infile;
	infile.open(metaDataPath.data());   //将文件流对象与文件连接起来 
	assert(infile.is_open());   //若失败,则输出错误消息,并终止程序运行 

	string s; //存储每一行的文本数据
	float ur_LAT = 0, ur_LON = 0, ll_LAR = 0, ll_LON = 0; //右上和左下经纬度
	while (getline(infile, s))
	{
		string S = removeNUllString(s);

		// 得到拍摄时间
		string sen = "DATE_ACQUIRED=";
		int time = S.find(sen);
		if (time != -1) {
			getTime = S.substr(sen.length(), S.length());
		}

		// 根据四个角度的经纬度，计算中央经纬度
		string UR_LAT = "CORNER_UR_LAT_PRODUCT=";
		time = S.find(UR_LAT);
		if (time != -1) {
			string temp = S.substr(UR_LAT.length(), S.length());
			ur_LAT = stringToNum<float>(temp);
		}

		string UR_LON = "CORNER_UR_LON_PRODUCT=";
		time = S.find(UR_LON);
		if (time != -1) {
			string temp = S.substr(UR_LON.length(), S.length());
			ur_LON = stringToNum<float>(temp);
		}

		string LL_LAR = "CORNER_LL_LAT_PRODUCT=";
		time = S.find(LL_LAR);
		if (time != -1) {
			string temp = S.substr(LL_LAR.length(), S.length());
			ll_LAR = stringToNum<float>(temp);
		}

		string LL_LON = "CORNER_LL_LON_PRODUCT=";
		time = S.find(LL_LON);
		if (time != -1) {
			string temp = S.substr(LL_LON.length(), S.length());
			ll_LON = stringToNum<float>(temp);
		}

		if ((ll_LAR != 0) && (ur_LAT != 0) && (ll_LON != 0) && (ur_LON != 0)) {
			*latitude = (ll_LAR + ur_LAT) / 2;
			*longtitude = (ll_LON + ur_LON) / 2;
			//printf("中央经度为:%f\n", *longtitude);
			//printf("中央纬度为:%f\n", *latitude);
		}

		// 得到 R波段辐射定标增益参数
		string senR = "RADIANCE_ADD_BAND_4=";
		time = S.find(senR);
		if (time != -1) {
			string temp = S.substr(senR.length(), S.length());
			*R_add = stringToNum<float>(temp);
			//printf("R波段辐射定标增益参数:%f\n",*R_add);
		}

		// 得到 NIR波段辐射定标增益参数
		string senNIR = "RADIANCE_ADD_BAND_5=";
		time = S.find(senNIR);
		if (time != -1) {
			string temp = S.substr(senNIR.length(), S.length());
			*NIR_add = stringToNum<float>(temp);
			//printf("NIR波段辐射定标增益参数:%f\n", *NIR_add);
		}

		// 得到 TIRS波段辐射定标增益参数
		string senTIRS = "RADIANCE_ADD_BAND_10=";
		time = S.find(senTIRS);
		if (time != -1) {
			string temp = S.substr(senTIRS.length(), S.length());
			*TIRS_add = stringToNum<float>(temp);
			//printf("TIRS波段辐射定标增益参数:%f\n", *TIRS_add);
		}

		// 得到 R波段辐射定标偏移参数
		string senRmult = "RADIANCE_MULT_BAND_4=";
		time = S.find(senRmult);
		if (time != -1) {
			string temp = S.substr(senRmult.length(), S.length());
			*R_mult = stringToNum<float>(temp);
			//printf("R波段辐射定标偏移参数:%f\n", *R_mult);
		}

		// 得到 NIR波段辐射定标偏移参数
		string senNIRmult = "RADIANCE_MULT_BAND_5=";
		time = S.find(senNIRmult);
		if (time != -1) {
			string temp = S.substr(senNIRmult.length(), S.length());
			*NIR_mult = stringToNum<float>(temp);
			//printf("NIR波段辐射定标偏移参数:%f\n", *NIR_mult);
		}

		// 得到 TIRS波段辐射定标偏移参数
		string senTIRSmult = "RADIANCE_MULT_BAND_10=";
		time = S.find(senTIRSmult);
		if (time != -1) {
			string temp = S.substr(senTIRSmult.length(), S.length());
			*TIRS_mult = stringToNum<float>(temp);
			//printf("TIRS波段辐射定标偏移参数:%f\n", *TIRS_mult);
		}

	}
	infile.close();
}

void Landsat8_LST::LandSat8Image_temperatureRetreval(char* RedBandPath, char* NirBandPath, char* TIRSBandPath,
	float t, float Lu, float Ld, string metaDataPath, char* imgPath) {

	// 读取红波段
	ImgData redband;
	readImage(RedBandPath, redband);

	// 读取近红波段
	ImgData nirband;
	readImage(NirBandPath, nirband);

	// 读取第10波段
	ImgData TIRSband;
	readImage(TIRSBandPath, TIRSband);

	// 读取影像元数据参数
	string getTime;
	float* R_add = new float;  //记得初始化局部变量指针
	float* NIR_add = new float;
	float* TIRS_add = new float;
	float* R_mult = new float;
	float* NIR_mult = new float;
	float* TIRS_mult = new float;
	float* longtitude = new float;
	float* latitude = new float;

	getMetaData(metaDataPath, R_add, NIR_add, TIRS_add, R_mult, NIR_mult,
		TIRS_mult, getTime, longtitude, latitude);

	// 辐射定标
	radImg((float*)(redband.pData), *R_add, *R_mult, redband.imglength); //R波段辐射定标
	radImg((float*)(nirband.pData), *NIR_add, *NIR_mult, redband.imglength); //NIR波段辐射定标
	radImg((float*)(TIRSband.pData), *TIRS_add, *TIRS_mult, redband.imglength); //热红外波段辐射定标

	// 比辐射率计算
	float* aima = SurEmiss(nirband.pData, redband.pData, redband.imglength);

	// 温度反演
	float* LandSurfaceTemp = LST(t, Lu, Ld, aima, TIRSband.pData, redband.imglength);

	// 导出结果
	writeImage(imgPath, LandSurfaceTemp, redband);

	delete[] redband.pData;
	delete[] nirband.pData;
	delete[] TIRSband.pData;
	delete R_add, NIR_add, TIRS_add, R_mult, NIR_mult, TIRS_mult, latitude;
}
